1 Introducción

En este trabajo se analiza la estructura multicapa de la red de transporte de la ciudad de Madrid a partir de dos conjuntos de datos reales proporcionados:

  • Un fichero con las paradas y estaciones de distintos modos de transporte.
  • Un fichero con información de viajes mensuales por distrito.

Los objetivos concretos son:

  • Construir una red multicapa donde cada capa representa un modo de transporte (bus, metro, bicimad, parking).
  • Estudiar la centralidad de las zonas de intercambio en cada capa.
  • Medir la versatilidad de las zonas (cómo cambia su importancia entre modos).
  • Realizar una visualización multicapa en la que se vean explícitamente las distintas capas conectadas.
  • Explorar de forma sencilla la vulnerabilidad del sistema ante fallos en una de las capas.

Toda la parte empírica se basa exclusivamente en los ficheros:

  • transporte_madrid_consolidado.csv
  • viajes_madrid.csv

2 Carga y exploración de datos

library(tidyverse)
library(janitor)
library(sf)
library(dbscan)
library(FNN)
library(igraph)
library(scales)
library(knitr)

2.1 Lectura de los ficheros

Se asume que los ficheros CSV están en el mismo directorio que este documento RMarkdown. Si no es así, ajustar las rutas.

ruta_transporte <- "transporte_madrid_consolidado.csv"
ruta_viajes     <- "viajes_madrid.csv"

transporte_raw <- read.csv(
  ruta_transporte,
  sep = ";",
)

viajes_raw <- read.csv(
  ruta_viajes,
  sep = ";",
)

glimpse(transporte_raw)
## Rows: 14,188
## Columns: 6
## $ stop_name         <chr> "                                                   …
## $ stop_lat          <dbl> 40.37835, 40.38341, 40.36813, 40.37773, 40.40753, 40…
## $ stop_lon          <dbl> -3.75458, -3.62481, -3.62283, -3.67311, -3.59176, -3…
## $ transport_mode    <chr> "bicimad                                            …
## $ stop_id           <chr> "382 -Parque Eugenia de Montijo, 2                  …
## $ distrito_asignado <int> 2807911, 2807918, 2807918, 2807913, 2807919, 2807911…
glimpse(viajes_raw)
## Rows: 8,064
## Columns: 11
## $ fecha             <chr> "2024-01", "2024-01", "2024-01", "2024-01", "2024-01…
## $ name              <chr> "Centro                                            "…
## $ edad              <chr> "0-25", "0-25", "0-25", "0-25", "0-25", "0-25", "0-2…
## $ sexo              <chr> "hombre", "hombre", "hombre", "hombre", "mujer", "mu…
## $ numero_viajes     <chr> "0", "1", "2", "2+", "0", "1", "2", "2+", "0", "1", …
## $ personas          <chr> "90541,756", "18809,207", "73248,991", "187551,534",…
## $ Codigo            <int> 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, …
## $ Provincia         <chr> "Madrid                ", "Madrid                ", …
## $ distrito          <dbl> 2807901, 2807901, 2807901, 2807901, 2807901, 2807901…
## $ zona_pernoctacion <int> 2807901, 2807901, 2807901, 2807901, 2807901, 2807901…
## $ poblacion         <int> 140644, 140644, 140644, 140644, 140644, 140644, 1406…

2.2 Descripción básica de transporte_madrid_consolidado

transporte <- clean_names(transporte_raw)

transporte$transport_mode <- str_trim(transporte$transport_mode)
transporte$transport_mode <- str_extract(transporte$transport_mode, "\\w+")

table(transporte$transport_mode)
## 
## bicimad     bus   metro parking 
##     626   12430    1050      82
columnas_interes <- transporte[, c("stop_lat", "stop_lon", "distrito_asignado")]
summary(columnas_interes)
##     stop_lat        stop_lon      distrito_asignado
##  Min.   :40.28   Min.   :-3.875   Min.   :2807901  
##  1st Qu.:40.40   1st Qu.:-3.711   1st Qu.:2807905  
##  Median :40.43   Median :-3.689   Median :2807909  
##  Mean   :40.43   Mean   :-3.683   Mean   :2807910  
##  3rd Qu.:40.46   3rd Qu.:-3.656   3rd Qu.:2807915  
##  Max.   :40.56   Max.   :-3.447   Max.   :2807921

Variables clave:

  • stop_name: nombre de la parada/estación.
  • stop_lat, stop_lon: coordenadas geográficas.
  • transport_mode: modo de transporte (bus, metro, bicimad, parking).
  • stop_id: identificador de parada.
  • distrito_asignado: código de distrito (INE).

Esto define los nodos base del sistema de transporte.

2.3 Descripción básica de viajes_madrid

library(readr)
library(dplyr)
library(janitor)
library(stringr)

viajes <- read.csv("viajes_madrid.csv",
                   sep = ";",
                   dec = ",")

viajes <- clean_names(viajes)

summary(viajes$personas)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4515   37486  110989  140469  186820  737342

Variables clave:

  • fecha: mes de referencia.
  • name: nombre del distrito.
  • edad, sexo: características sociodemográficas.
  • numero_viajes: categoría de nº de viajes (0, 1, 2, 2+).
  • personas: número de personas en esa celda.
  • distrito: código de distrito.
  • poblacion: población total del distrito.

2.3.1 Indicador de intensidad de viajes por distrito

Construimos un indicador sencillo de viajes medios por persona en cada distrito. Primero pasamos la categoría numero_viajes a un valor numérico aproximado y luego calculamos una media ponderada por el número de personas.

2.3.2 Indicador de viajes medios diarios por persona y distrito

La variable personas se encontraba codificada con formato decimal europeo, por lo que fue necesario convertirla explícitamente a formato numérico. Asimismo, la variable numero_viajes es categórica (0, 1, 2, 2+), por lo que se transformó en una aproximación numérica asignando el valor 3 a la categoría “2+”. El indicador final de viajes medios por persona se calculó como una media ponderada por el número de personas en cada categoría, obteniéndose valores en torno a 2 viajes diarios por persona, coherentes con patrones conocidos de movilidad urbana.

library(dplyr)
library(knitr)

viajes$viajes_cat <- case_when(
  viajes$numero_viajes == "0"  ~ 0,
  viajes$numero_viajes == "1"  ~ 1,
  viajes$numero_viajes == "2"  ~ 2,
  viajes$numero_viajes == "2+" ~ 3,
  TRUE ~ NA_real_
)

viajes_distrito <- aggregate(
  cbind(personas, viajes_totales = viajes_cat * personas) ~ distrito, 
  data = viajes, 
  FUN = sum, 
  na.rm = TRUE
)

viajes_distrito$viajes_medios_por_persona <- viajes_distrito$viajes_totales / viajes_distrito$personas
viajes_distrito$personas_totales <- viajes_distrito$personas

kable(
  head(viajes_distrito),
  digits = 3,
  caption = "Viajes medios diarios por persona y distrito"
)
Viajes medios diarios por persona y distrito
distrito personas viajes_totales viajes_medios_por_persona personas_totales
2807901 51888948 110673624 2.133 51888948
2807902 51542198 108276834 2.101 51542198
2807903 39251902 78164164 1.991 39251902
2807904 49039609 98112207 2.001 49039609
2807905 48349024 95741527 1.980 48349024
2807906 53679153 110757611 2.063 53679153
summary(viajes_distrito$viajes_medios_por_persona)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.956   2.004   2.038   2.041   2.063   2.135

Interpretación de los Resultados (Estadísticas de Viajes)

El resumen estadístico (summary) de la variable viajes_medios_por_persona revela un patrón de comportamiento extremadamente homogéneo en todos los distritos analizados:

  1. Uniformidad del Comportamiento: El rango de valores es muy estrecho, oscilando apenas entre un mínimo de 1.956 y un máximo de 2.135 viajes por persona. Esto indica que, independientemente del distrito (sea centro o periferia), la movilidad media individual es muy constante.

  2. Tendencia Central (La regla del “Ida y Vuelta”): Tanto la media (2.041) como la mediana (2.038) se sitúan prácticamente en 2. Esto tiene un sentido lógico muy fuerte en transporte urbano: representa el patrón estándar de un viaje de ida y otro de vuelta (péndulo diario).

  3. Implicación para el Modelo: Dado que la desviación entre el primer cuartil (2.004) y el tercer cuartil (2.063) es mínima, esta variable no actúa como un factor discriminante fuerte por sí sola. Esto significa que las diferencias de flujo total que observaremos en la red no se deben a que la gente de un barrio viaje “más veces” que la de otro, sino probablemente a la densidad de población de cada zona (mismo ratio de viajes por persona, pero mucha más gente en unas zonas que en otras).

Se observa una notable homogeneidad en la movilidad individual a lo largo de Madrid, con una media de 2.04 viajes por persona y una varianza despreciable. Esto sugiere que la intensidad de uso de la red de transporte depende casi exclusivamente del volumen demográfico de cada distrito y no de diferencias significativas en los hábitos de movilidad individual entre zonas.

3 Construcción de la red multicapa

La idea es pasar de paradas individuales a zonas de intercambio (clusters de paradas cercanas), y luego construir una capa por modo de transporte.

3.1 1. Definición de zonas de intercambio (clustering espacial)

# Creamos objeto sf con coordenadas
transporte_sf <- st_as_sf(transporte, coords = c("stop_lon", "stop_lat"), crs = 4326, remove = FALSE)
transporte_sf <- st_transform(transporte_sf, 25830)  # sistema métrico (UTM zona 30)

coords_m <- st_coordinates(transporte_sf)
coords_m <- as.matrix(coords_m)

# Clustering DBSCAN: paradas a menos de 150 m se agrupan
set.seed(123)
db <- dbscan(coords_m, eps = 150, minPts = 2)

transporte_sf$cluster_id <- db$cluster

# Tratamos el ruido (cluster 0) como clusters individuales
max_id <- max(transporte_sf$cluster_id)
es_ruido <- transporte_sf$cluster_id == 0
transporte_sf$cluster_id[es_ruido] <- max_id + seq_len(sum(es_ruido))

# Zona de intercambio = centroide del cluster
df_no_geo <- st_drop_geometry(transporte_sf)

# Calculamos medias de coordenadas
zonas <- aggregate(cbind(stop_lat, stop_lon) ~ cluster_id, 
                   data = df_no_geo, 
                   FUN = mean)

# distrito asignado más frecuente en el cluster
distritos_freq <- aggregate(distrito_asignado ~ cluster_id, 
                             data = df_no_geo, 
                             FUN = function(x) as.integer(names(which.max(table(x)))))

# Modos disponibles
modos_agg <- aggregate(transport_mode ~ cluster_id, 
                       data = df_no_geo, 
                       FUN = function(x) paste(sort(unique(x)), collapse = ","))

# Unimos todo en el objeto zonas
zonas <- merge(zonas, distritos_freq, by = "cluster_id")
zonas <- merge(zonas, modos_agg, by = "cluster_id")

nrow(zonas)
## [1] 1558
head(zonas)

Cada cluster_id representa una zona de intercambio multimodal.

Representamos la distribución espacial aproximada (proyección geográfica simple):

# Aseguramos que sea factor
zonas_sf <- st_as_sf(zonas, coords = c("stop_lon", "stop_lat"), crs = 4326)

zonas_sf$distrito_asignado <- as.factor(zonas_sf$distrito_asignado)
n_distritos <- length(levels(zonas_sf$distrito_asignado))

# Usamos hcl.colors con la paleta "Dark 3" o "Set 2", que son mucho más elegantes
plot(zonas_sf["distrito_asignado"], 
     main = "Zonas de intercambio por distrito",
     pal = hcl.colors(n_distritos, palette = "Dark 3"), 
     key.pos = 4,
     pch = 19, 
     cex = 0.6,
     border = NA) # Quitamos bordes si los hubiera para que el color sea más puro

3.2 2. Nodos por modo de transporte

modos <- sort(unique(transporte_sf$transport_mode))

# Creamos una tabla de frecuencias y la convertimos a data frame
df_tab <- as.data.frame.matrix(table(transporte_sf$cluster_id, transporte_sf$transport_mode))

# Convertimos las frecuencias en presencia (1) o ausencia (0)
df_tab[df_tab > 0] <- 1
df_tab$cluster_id <- as.numeric(rownames(df_tab))

# Unimos con las zonas
zonas_attrs <- merge(zonas, df_tab, by = "cluster_id", all.x = TRUE)

# Reemplazamos NAs por 0 en las columnas de modos
zonas_attrs[modos][is.na(zonas_attrs[modos])] <- 0

kable(head(zonas_attrs), caption = "Atributos básicos de las zonas por modo")
Atributos básicos de las zonas por modo
cluster_id stop_lat stop_lon distrito_asignado transport_mode bicimad bus metro parking
1 40.38368 -3.624114 2807918 bicimad,bus,metro 1 1 1 0
2 40.38096 -3.728060 2807911 bicimad,bus,metro 1 1 1 0
3 40.38626 -3.719795 2807912 bicimad,bus,metro 1 1 1 0
4 40.37886 -3.732537 2807911 bicimad,bus 1 1 0 0
5 40.40855 -3.672909 2807903 bicimad,bus 1 1 0 0
6 40.40380 -3.706547 2807901 bicimad,metro 1 0 1 0

En este análisis, la estructura de los clústeres utiliza la proximidad geográfica para configurar nodos de naturaleza multimodal. Cada clúster integra diversas tipologías de infraestructuras de transporte basándose en un criterio de agrupación por distancias mínimas, lo que permite consolidar puntos de acceso heterogéneos en una unidad funcional única.

Como se observa en la matriz resultante, la presencia de servicios como \(bicimad\), \(bus\), \(metro\) y \(parking\) dentro de un mismo \(cluster\_id\) evidencia que el modelo no solo identifica cercanía espacial, sino que captura la interconectividad del tejido urbano, facilitando la transición entre diferentes modos de movilidad en puntos estratégicos de la red.

3.3 3. Construcción de redes intra-capa

Este bloque de código construye la estructura topológica intra-capa (conexiones internas) para cada medio de transporte mediante un criterio de proximidad espacial. En primer lugar, proyecta las coordenadas geográficas a un sistema métrico (UTM) para poder calcular distancias reales en metros. Posteriormente, aplica un algoritmo de K-Vecinos Más Cercanos (KNN) que conecta cada parada (nodo) con sus 3 estaciones más próximas del mismo modo de transporte, asignando a cada conexión un peso (\(w\)) inversamente proporcional a la distancia física; de esta forma, se generan redes de vecindad geométrica donde los nodos cercanos tienen enlaces más fuertes, simulando la accesibilidad local dentro de cada capa.

# 1. Preparación de coordenadas de forma sencilla
centroides_mat <- st_coordinates(st_transform(zonas_sf, 25830))
zonas_centroides <- data.frame(
  cluster_id = zonas$cluster_id,
  x = centroides_mat[, 1],
  y = centroides_mat[, 2]
)

# 2. Función con lógica explícita
build_edges_simple <- function(modo_nombre, k = 3) {
  
  # Filtrar qué zonas tienen este modo activo
  ids_con_modo <- zonas_attrs$cluster_id[zonas_attrs[[modo_nombre]] == 1]
  
  # Si no hay suficientes zonas para conectar, devolvemos tabla vacía
  if (length(ids_con_modo) < 2) return(data.frame())
  
  # Obtener coordenadas solo de esas zonas
  coords <- zonas_centroides[zonas_centroides$cluster_id %in% ids_con_modo, ]
  
  lista_aristas <- list()
  
  # Bucle: Para cada zona, calculamos la distancia a TODAS las demás
  for (i in 1:nrow(coords)) {
    punto_origen <- coords[i, ]
    
    # Calculamos distancia euclidiana: sqrt((x1-x2)^2 + (y1-y2)^2)
    distancias <- sqrt((punto_origen$x - coords$x)^2 + (punto_origen$y - coords$y)^2)
    
    # Creamos un dataframe temporal con las distancias de 'i' al resto
    df_dist <- data.frame(
      modo = modo_nombre,
      from = punto_origen$cluster_id,
      to   = coords$cluster_id,
      dist = distancias
    )
    
    # Quitamos la distancia a sí mismo (que es 0)
    df_dist <- df_dist[df_dist$from != df_dist$to, ]
    
    # Ordenamos por distancia y nos quedamos con los 'k' más cercanos
    df_dist <- df_dist[order(df_dist$dist), ]
    df_dist <- head(df_dist, k)
    
    lista_aristas[[i]] <- df_dist
  }
  
  # Unir todos los resultados del bucle en una tabla
  resultado <- do.call(rbind, lista_aristas)
  
  # Calcular el peso (w) de forma manual
  resultado$w <- 1 / (1 + resultado$dist)
  
  return(resultado)
}

edges_intra_lista <- list()
for (m in modos) {
  edges_intra_lista[[m]] <- build_edges_simple(m)
}

# Combinamos todas las aristas en un solo dataframe
edges_intra <- do.call(rbind, edges_intra_lista)

kable(head(edges_intra), caption = "Ejemplo de aristas intra-capa")
Ejemplo de aristas intra-capa
modo from to dist w
bicimad.144 bicimad 1 288 649.2619 0.0015378
bicimad.384 bicimad 1 1382 712.9355 0.0014007
bicimad.218 bicimad 1 515 746.3214 0.0013381
bicimad.149 bicimad 2 300 321.5842 0.0031000
bicimad.4 bicimad 2 4 445.7282 0.0022385
bicimad.428 bicimad 2 1463 736.8067 0.0013554

Descripción del Dataframe obtenido (edges_intra)

El objeto resultante es una tabla consolidada (tibble) donde cada fila representa una arista (conexión) entre dos paradas. Sus columnas son:

  • modo: Identificador de la capa de transporte (ej. “metro”, “bus”, “bicimad”).

  • from: ID numérico del nodo de origen (cluster_id).to: ID numérico del nodo de destino (uno de sus vecinos más cercanos).

  • dist: La distancia euclídea en metros entre ambos nodos.

  • w: El peso de la conexión, calculado como \(1/(1 + \text{distancia})\). Varía entre 0 y 1, siendo mayor cuanto más cerca están las paradas.

3.4 4. Grafos por capa (lista g_list)

Este bloque tiene una misión crítica en el análisis multicapa: garantizar la alineación de nodos (Node Alignment). La función build_graph_modo transforma la lista de enlaces (edges_intra) en objetos de grafo (igraph), pero con una particularidad vital: asegura que todas las capas tengan exactamente los mismos nodos, incluso si en una capa específica un nodo no tiene conexiones (está aislado). El código identifica los nodos que faltan en una capa (setdiff) y los añade manualmente como vértices aislados. Sin este paso, las matrices de adyacencia tendrían dimensiones distintas y sería matemáticamente imposible construir el objeto “multiplex” (la “lasaña” no encajaría).

nodos_ids <- zonas_attrs$cluster_id

build_graph_simple <- function(nombre_modo) {
  filas_modo <- edges_intra[edges_intra$modo == nombre_modo, ]
  e_modo <- filas_modo[, c("from", "to", "w")]

  g <- graph_from_data_frame(
    d = e_modo, 
    directed = FALSE, 
    vertices = zonas_attrs
  )

  nodos_en_grafo <- as.integer(V(g)$name)
  faltantes <- nodos_ids[!(nodos_ids %in% nodos_en_grafo)]
  
  if (length(faltantes) > 0) {
    g <- g + vertices(as.character(faltantes))
  }

  return(g)
}

g_list <- list()
for (m in modos) {
  g_list[[m]] <- build_graph_simple(m)
}

g_list
## $bicimad
## IGRAPH 626af4c UN-- 1558 1404 -- 
## + attr: name (v/c), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), transport_mode (v/c), bicimad (v/n), bus (v/n), metro (v/n),
## | parking (v/n), w (e/n)
## + edges from 626af4c (vertex names):
##  [1] 1 --288  1 --1382 1 --515  2 --300  2 --4    2 --1463 3 --1511 3 --1241
##  [9] 3 --1547 2 --4    4 --300  4 --35   5 --283  5 --1108 5 --240  6 --927 
## [17] 6 --1496 6 --239  7 --29   7 --479  7 --239  8 --735  8 --1342 8 --756 
## [25] 11--248  11--41   11--460  12--372  12--488  12--1070 13--1391 13--112 
## [33] 13--436  14--1556 14--1056 14--156  16--583  16--1350 16--728  23--778 
## [41] 23--1354 23--1398 24--626  24--306  24--248  25--296  25--748  25--370 
## + ... omitted several edges
## 
## $bus
## IGRAPH 626ba14 UN-- 1558 4038 -- 
## + attr: name (v/c), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), transport_mode (v/c), bicimad (v/n), bus (v/n), metro (v/n),
## | parking (v/n), w (e/n)
## + edges from 626ba14 (vertex names):
##  [1] 1 --1289 1 --280  1 --1231 2 --1026 2 --300  2 --1304 3 --1303 3 --1082
##  [9] 3 --1148 4 --652  4 --532  4 --229  5 --283  5 --1108 5 --945  7 --29  
## [17] 7 --479  7 --239  8 --819  8 --830  8 --1153 9 --390  9 --335  9 --48  
## [25] 10--844  10--354  10--1446 11--767  11--248  11--41   12--1288 12--368 
## [33] 12--1449 13--638  13--639  13--642  14--158  14--1272 14--1198 15--422 
## [41] 15--357  15--187  16--969  16--726  16--724  18--20   18--19   18--345 
## + ... omitted several edges
## 
## $metro
## IGRAPH 626bf74 UN-- 1558 684 -- 
## + attr: name (v/c), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), transport_mode (v/c), bicimad (v/n), bus (v/n), metro (v/n),
## | parking (v/n), w (e/n)
## + edges from 626bf74 (vertex names):
##  [1] 1 --230  1 --227  1 --49   2 --229  2 --1077 2 --3    2 --3    3 --459 
##  [9] 3 --816  6 --239  6 --7    6 --132  7 --239  6 --7    7 --118  11--248 
## [17] 11--41   11--24   12--372  12--746  12--510  13--112  13--436  13--97  
## [25] 17--187  17--18   17--251  18--187  18--30   18--68   24--248  24--838 
## [33] 11--24   27--917  27--572  27--116  28--287  28--696  28--523  30--68  
## [41] 18--30   30--58   40--145  40--164  40--245  41--753  11--41   40--41  
## + ... omitted several edges
## 
## $parking
## IGRAPH 626c4d2 UN-- 1558 162 -- 
## + attr: name (v/c), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), transport_mode (v/c), bicimad (v/n), bus (v/n), metro (v/n),
## | parking (v/n), w (e/n)
## + edges from 626c4d2 (vertex names):
##  [1] 7  --29   7  --1488 7  --239  8  --1532 8  --41   8  --259  27 --766 
##  [8] 27 --815  27 --116  28 --1475 28 --292  28 --815  29 --239  7  --29  
## [15] 29 --657  41 --245  41 --1503 41 --164  55 --1052 55 --657  55 --423 
## [22] 56 --814  56 --413  56 --67   67 --814  56 --67   67 --121  71 --446 
## [29] 71 --1052 55 --71   84 --657  84 --1052 84 --239  101--164  101--573 
## [36] 101--102  102--573  102--413  102--446  116--766  116--815  27 --116 
## + ... omitted several edges

La ejecución del algoritmo de construcción de grafos ha resultado en la generación de cuatro objetos igraph independientes, correspondientes a las capas de BiciMAD, Bus, Metro y Parking. La inspección de la salida (g_list) revela las siguientes características estructurales fundamentales para el modelo multicapa:

  1. Isomorfismo Nodal (Alineación de Capas) Se confirma que las cuatro capas poseen exactamente 1558 nodos (vertices = 1558). Este resultado valida el procedimiento de alineación de nodos (node alignment), garantizando que todas las capas comparten el mismo conjunto de definiciones espaciales (paradas/zonas), incluso si un modo de transporte no tiene actividad en una zona específica. Esto es un requisito matemático indispensable para la construcción del tensor supra-adyacente en muxViz.

  2. Heterogeneidad en la Densidad de Conexiones Se observa una gran disparidad en el número de aristas (edges) entre capas, lo que refleja la naturaleza operativa de cada modo:

  • Bus (4038 enlaces): Es la capa con mayor densidad topológica. Su elevado número de aristas indica una alta capilaridad, actuando como el tejido conectivo principal que une distancias cortas y medias con una gran redundancia local.

  • BiciMAD (1404 enlaces): Presenta una conectividad intermedia. Aunque inferior al autobús, su densidad sugiere una fuerte cohesión en las zonas donde está implantado (principalmente la almendra central), formando clústeres locales densos.

  • Metro (684 enlaces): Muestra una topología más dispersa. Al funcionar como una red “esqueleto” o vertebral, sus conexiones son menores en número pero estratégicas, diseñadas para el transporte rápido de larga distancia sin saturar el espacio con conexiones redundantes.

  • Parking (162 enlaces): Es la capa más desconectada. La escasez de aristas refleja que los aparcamientos funcionan operativamente como “nodos satélite” o destinos finales, con muy poca interacción directa entre ellos bajo el criterio de vecindad espacial establecido (\(k=3\)).

  1. Preservación de Atributos La salida confirma (attr: name, stop_lat…) que se han preservado los metadatos espaciales y demográficos en todos los nodos. Esto permitirá que, en las fases posteriores de visualización y análisis de vulnerabilidad, se puedan correlacionar las métricas de red (como el PageRank) con variables del mundo real como el distrito asignado o los viajes totales.
Tabla 1: Resumen Topológico de las Capas de la Red Multicapa
Capa Nodos Aristas Densidad_Relativa Funcion_Sistema
Bus 1558 4038 Alta Capilaridad / Cobertura
BiciMAD 1558 1404 Media Movilidad última milla / Centro
Metro 1558 684 Baja Vertebración / Eficiencia
Parking 1558 162 Muy Baja Nodos Satélite / Intercambio

4 Visualización multicapa con muxViz

En esta sección generamos la visualización clave: un gráfico 3D en el que se ven las capas (modos) apiladas, con los nodos alineados entre capas.

Para ello usamos la función plot_multiplex3D de muxViz. Se asume que muxViz está instalado; si no lo está, el código no fallará, pero mostrará un mensaje.

# --- 1. PREPARACIÓN DE DATOS ---
zonas_attrs <- zonas_attrs[order(zonas_attrs$cluster_id), ]

centroides_mat <- st_coordinates(st_transform(zonas_sf, 25830))
zonas_centroides <- data.frame(
  cluster_id = zonas$cluster_id,
  x = centroides_mat[, 1],
  y = centroides_mat[, 2]
)
zonas_centroides <- zonas_centroides[order(zonas_centroides$cluster_id), ]


# --- 2. GENERACIÓN DE ARISTAS (INTRA-CAPA) ---
edges_intra_lista <- list()

for (m in modos) {
  ids_modo <- zonas_attrs$cluster_id[zonas_attrs[[m]] == 1]
  
  if (length(ids_modo) >= 2) {
    coords_sub <- zonas_centroides[zonas_centroides$cluster_id %in% ids_modo, ]
    aristas_modo <- list()
    
    for (i in 1:nrow(coords_sub)) {
      dists <- sqrt((coords_sub$x[i] - coords_sub$x)^2 + (coords_sub$y[i] - coords_sub$y)^2)
      
      df_temp <- data.frame(
        modo = m,
        from = coords_sub$cluster_id[i],
        to   = coords_sub$cluster_id,
        dist = dists
      )
      
      df_temp <- df_temp[df_temp$from != df_temp$to, ]
      df_temp <- df_temp[order(df_temp$dist), ]
      aristas_modo[[i]] <- head(df_temp, 3)
    }
    edges_intra_lista[[m]] <- do.call(rbind, aristas_modo)
  }
}

edges_intra <- do.call(rbind, edges_intra_lista)
edges_intra$w <- 1 / (1 + edges_intra$dist)


# --- 3. CONSTRUCCIÓN DE GRAFOS ---
g_list <- list()

for (m in modos) {
  e_sub <- edges_intra[edges_intra$modo == m, c("from", "to", "w")]
  g <- graph_from_data_frame(d = e_sub, directed = FALSE, vertices = zonas_attrs)
  g_list[[m]] <- g
}


# --- 4. VISUALIZACIÓN 3D MEJORADA (MÁS ESPACIO) ---
if (requireNamespace("rgl", quietly = TRUE)) {
  library(rgl)
  library(igraph)
  
  clear3d() 
  open3d()
  par3d(windowRect = c(0, 0, 900, 900)) # Ventana un poco más grande
  bg3d("white")
  view3d(theta = 30, phi = 30, zoom = 0.7) # Zoom un poco más lejos
  
  # --- AJUSTES DE ESPACIADO ---
  # 1. Factor de dispersión horizontal: Multiplicamos para separar los puntos XY
  spread_factor <- 3.0 
  lay_norm <- scale(centroides_mat) * spread_factor
  
  # 2. Separación vertical aumentada
  z_sep <- 8.0 
  
  colores <- c(bicimad="#1a9641", bus="#d7191c", metro="#2c7bb6", parking="#fdae61")
  layer_names <- names(g_list)
  
  # --- BUCLE: CAPA POR CAPA ---
  for (i in 1:length(g_list)) {
    g <- g_list[[i]]
    lay_z <- rep(i * z_sep, nrow(lay_norm)) 
    
    # A. PLANO BASE
    rango_x <- range(lay_norm[,1])
    rango_y <- range(lay_norm[,2])
    # Hacemos el plano ligeramente más grande que los puntos
    quads3d(
      x = c(rango_x[1]-1, rango_x[2]+1, rango_x[2]+1, rango_x[1]-1),
      y = c(rango_y[1]-1, rango_y[1]-1, rango_y[2]+1, rango_y[2]+1),
      z = rep(i * z_sep - 0.2, 4), 
      col = "gray95", alpha = 0.1 # Más transparente aún
    )
    
    # B. ARISTAS
    edgelist <- as_edgelist(g, names = FALSE)
    if (nrow(edgelist) > 0) {
      x0 <- lay_norm[edgelist[,1], 1]; x1 <- lay_norm[edgelist[,2], 1]
      y0 <- lay_norm[edgelist[,1], 2]; y1 <- lay_norm[edgelist[,2], 2]
      z0 <- rep(i * z_sep, length(x0))
      
      segments3d(
        x = as.vector(rbind(x0, x1)),
        y = as.vector(rbind(y0, y1)),
        z = as.vector(rbind(z0, z0)),
        col = "gray0", lwd = 1, alpha = 0.4
      )
    }
    
    # C. NODOS
    deg <- degree(g)
    active_idx <- which(deg > 0)
    
    if (length(active_idx) > 0) {
      # Ajuste fino del tamaño de las esferas
      sizes <- 0.8 + (deg[active_idx] / (max(deg)+1)) * 1.5
      spheres3d(
        x = lay_norm[active_idx, 1], 
        y = lay_norm[active_idx, 2], 
        z = lay_z[active_idx],
        radius = sizes * 0.2, # Radio ajustado
        col = colores[i]
      )
    }
    
    # D. ETIQUETA
    text3d(
      x = min(lay_norm[,1]) - 1, # Un poco más a la izquierda
      y = min(lay_norm[,2]) - 1, 
      z = i * z_sep + 0.5, 
      text = layer_names[i], col = colores[i], cex = 1.2, font = 2
    )
  }
  
  # --- CONEXIONES VERTICALES ---
  for (k in 1:nrow(lay_norm)) {
    # Comprobamos grado en cada grafo de la lista
    capas_activas <- sapply(g_list, function(gr) degree(gr)[k] > 0)
    idx_capas <- which(capas_activas)
    
    if (length(idx_capas) >= 2) {
      z_vals <- idx_capas * z_sep
      segments3d(
        x = rep(lay_norm[k,1], 2),
        y = rep(lay_norm[k,2], 2),
        z = range(z_vals),
        col = "gray0", lwd = 1, alpha = 0.15 # Más sutiles para no ensuciar
      )
    }
  }
  
  rglwidget()
  
} else {
  message("Instala rgl")
}

Visualización de la red por capas y agregada

# --- CONFIGURACIÓN ESTÉTICA ---
colores <- c(bicimad="#1a9641", bus="#d7191c", metro="#2c7bb6", parking="#fdae61")
titulos <- names(g_list)
# Layout geográfico fijo (usamos las coordenadas que ya calculaste)
layout_geo <- as.matrix(centroides_mat)


# ==============================================================================
# VISTA 1: PANEL 2D (MULTIPLEX SLICES)
# Muestra cada capa por separado manteniendo la posición geográfica
# ==============================================================================
par(mfrow = c(2, 2), mar = c(1, 1, 3, 1)) # Rejilla 2x2

for (m in titulos) {
  g <- g_list[[m]]
  deg <- degree(g)
  
  # Lógica visual: 
  # Nodos activos = color intenso y tamaño según grado
  # Nodos inactivos = gris muy claro y pequeños (para mantener referencia espacial)
  v_col <- ifelse(deg > 0, colores[m], adjustcolor("gray90", alpha.f = 0.5))
  v_size <- ifelse(deg > 0, 3 + (deg/max(deg))*5, 1)
  v_frame <- ifelse(deg > 0, "black", NA) # Sin borde si está inactivo
  
  # Aristas: Más finas y transparentes
  e_col <- adjustcolor(colores[m], alpha.f = 0.6)
  
  plot(g,
       layout = layout_geo,
       vertex.label = NA,
       vertex.size = v_size,
       vertex.color = v_col,
       vertex.frame.color = v_frame,
       edge.color = e_col,
       edge.width = 0.8,
       edge.arrow.size = 0,
       main = toupper(m))
  
  box(col = "gray80") # Marco suave alrededor de cada plot
}

# ==============================================================================
# VISTA 2: RED AGREGADA (OVERLAY)
# Superpone todo para ver "corredores" multimodales
# ==============================================================================
par(mfrow = c(1, 1), mar = c(0, 0, 2, 0))

# Creamos un plot vacío primero para establecer los límites
plot(layout_geo, type = "n", axes = FALSE, xlab = "", ylab = "", 
     main = "RED AGREGADA (OVERLAP)")

# Dibujamos las aristas de cada capa acumulativamente
# Truco: Dibujamos primero las capas más densas (Bus/Metro) al fondo, Bicimad arriba
orden_capas <- c("bus", "metro", "parking", "bicimad") 

for (m in orden_capas) {
  if (m %in% names(g_list)) {
    g <- g_list[[m]]
    
    # Extraemos aristas solo si existen
    if (gsize(g) > 0) {
      el <- as_edgelist(g, names = FALSE)
      
      # Coordenadas de inicio y fin de las aristas
      x0 <- layout_geo[el[,1], 1]
      y0 <- layout_geo[el[,1], 2]
      x1 <- layout_geo[el[,2], 1]
      y1 <- layout_geo[el[,2], 2]
      
      # Dibujar líneas (curvadas ligeramente para ver solapamientos sería ideal, 
      # pero rectas es más limpio en mapas densos)
      segments(x0, y0, x1, y1, 
               col = adjustcolor(colores[m], alpha.f = 0.4), 
               lwd = 1)
    }
  }
}

# Dibujamos nodos ENCIMA de las líneas
# Calculamos grado total sumando todas las capas
deg_matrix <- sapply(g_list, degree)
deg_total <- rowSums(deg_matrix)
active_nodes <- deg_total > 0

points(layout_geo[active_nodes, ], 
       pch = 21, 
       bg = "white", 
       col = "gray30", 
       cex = 0.8 + (deg_total[active_nodes]/max(deg_total))*1.5)

legend("topright", legend = names(colores), 
       col = colores, lty = 1, lwd = 3, bty = "n", cex = 0.8)

Finalmente, para poder tener una perspectiva distinta utilizamos el resto de la visualización en un entorno tridimensional interactivo. A diferencia de las proyecciones planas, este gráfico ‘desapila’ la red de transporte ubicando cada modo (Metro, Bus, BiciMAD, Parking) en un nivel diferente del eje Z. Esto nos permite observar la estructura multiplex del sistema: las líneas verticales negras conectan una misma zona geográfica a través de sus diferentes capas, facilitando la identificación visual de los nodos con mayor integración multimodal.

# ===========================================
# Visualización 3D multicapa (inter-capa con scatter3d lines)
# ===========================================

if (requireNamespace("plotly", quietly = TRUE)) {

  library(plotly)
  library(scales)

  # ---- 1. Capas y eje Z
  layer_order <- modos
  layer_index <- setNames(seq_along(layer_order) - 1, layer_order)  # 0,1,2,3

  # ---- 2. Nodos multilayer: solo zonas donde ese modo existe
  nodes_list <- list()
  for (m in layer_order) {
    # Filtramos zonas con ese modo activo
    sub_zonas <- zonas_attrs[zonas_attrs[[m]] == 1, ]
    
    # Cruzamos con centroides
    df_m <- merge(sub_zonas[, "cluster_id", drop = FALSE], zonas_centroides, by = "cluster_id")
    
    # Construimos el dataframe de la capa
    capa_df <- data.frame(
      node       = paste(df_m$cluster_id, m, sep = "_"),
      cluster_id = df_m$cluster_id,
      layer      = m,
      x_raw      = df_m$x,
      y_raw      = df_m$y,
      stringsAsFactors = FALSE
    )
    nodes_list[[m]] <- capa_df
  }
  nodes_multi_3d <- do.call(rbind, nodes_list)

  # Normalizamos X, Y y asignamos Z
  nodes_multi_3d$x <- rescale(nodes_multi_3d$x_raw, to = c(0, 1))
  nodes_multi_3d$y <- rescale(nodes_multi_3d$y_raw, to = c(0, 1))
  nodes_multi_3d$z <- layer_index[nodes_multi_3d$layer]

  # ---- 3. Multimodalidad por zona
  # Calculamos cuántos modos hay activos sumando las columnas de modos
  modos_activos_vec <- rowSums(zonas_attrs[, modos])
  node_multimodalidad <- data.frame(
    cluster_id = zonas_attrs$cluster_id,
    modos_activos = modos_activos_vec
  )

  # Clusters con al menos 2 modos
  clusters_interes <- node_multimodalidad$cluster_id[node_multimodalidad$modos_activos >= 2]

  # ---- 4. Aristas inter-capa: todos los pares de modos presentes
  edges_inter_list <- list()
  for (cid in clusters_interes) {
    fila <- zonas_attrs[zonas_attrs$cluster_id == cid, modos]
    modos_presentes <- modos[as.logical(unlist(fila))]

    if (length(modos_presentes) >= 2) {
      comb <- t(combn(modos_presentes, 2))
      edges_inter_list[[as.character(cid)]] <- data.frame(
        from       = paste(cid, comb[, 1], sep = "_"),
        to         = paste(cid, comb[, 2], sep = "_"),
        cluster_id = cid,
        stringsAsFactors = FALSE
      )
    }
  }
  edges_inter_3d <- do.call(rbind, edges_inter_list)

  cat("Número de aristas inter-capa generadas:", nrow(edges_inter_3d), "\n")

  # ---- 5. Posiciones de nodos para las aristas
  pos_nodes <- nodes_multi_3d[, c("node", "x", "y", "z", "layer")]

  edges_plot <- merge(edges_inter_3d, pos_nodes, by.x = "from", by.y = "node")
  edges_plot <- merge(edges_plot, pos_nodes, by.x = "to", by.y = "node", suffixes = c("", "_to"))
  
  # Renombrar manualmente para claridad
  names(edges_plot)[names(edges_plot) == "x_to"] <- "xend"
  names(edges_plot)[names(edges_plot) == "y_to"] <- "yend"
  names(edges_plot)[names(edges_plot) == "z_to"] <- "zend"

  cat("Aristas con coordenadas válidas:", nrow(edges_plot), "\n")

  # ---- 6. Convertimos aristas a formato "long" para scatter3d lines
  if (nrow(edges_plot) > 0) {
    # Creamos vectores con NAs intermedios para separar líneas en plotly
    x_long <- as.vector(t(cbind(edges_plot$x, edges_plot$xend, NA)))
    y_long <- as.vector(t(cbind(edges_plot$y, edges_plot$yend, NA)))
    z_long <- as.vector(t(cbind(edges_plot$z, edges_plot$zend, NA)))
    edges_long <- data.frame(x = x_long, y = y_long, z = z_long)
  } else {
    edges_long <- data.frame(x = numeric(), y = numeric(), z = numeric())
  }

  cat("Puntos totales en edges_long:", nrow(edges_long), "\n")

  # ---- 7. Añadimos multimodalidad y tamaño de nodo
  nodes_multi_3d <- merge(nodes_multi_3d, node_multimodalidad, by = "cluster_id", all.x = TRUE)
  nodes_multi_3d$modos_activos[is.na(nodes_multi_3d$modos_activos)] <- 1
  nodes_multi_3d$size <- rescale(nodes_multi_3d$modos_activos, to = c(4, 11))

  # ---- 8. Colores por capa
  layer_colors <- c(
    bicimad = "#1a9641",
    bus     = "#d7191c",
    metro   = "#2c7bb6",
    parking = "#fdae61"
  )

  # ---- 9. Gráfico 3D
  p3d <- plot_ly()

  # Nodos
  p3d <- add_trace(
    p3d,
    data = nodes_multi_3d,
    type = "scatter3d",
    mode = "markers",
    x = ~x, y = ~y, z = ~z,
    color = ~layer,
    colors = layer_colors,
    marker = list(size = ~size),
    hoverinfo = "text",
    text = ~paste(
      "Zona:", cluster_id,
      "<br>Modo:", layer,
      "<br>Modos activos (total cluster):", modos_activos
    )
  )

  # Aristas inter-capa como líneas negras
  if (nrow(edges_long) > 0) {
    p3d <- add_trace(
      p3d,
      data = edges_long,
      type = "scatter3d",
      mode = "lines",
      x = ~x, y = ~y, z = ~z,
      line = list(color = "black", width = 4),
      showlegend = FALSE,
      hoverinfo = "none"
    )
  }

  p3d <- layout(
    p3d,
    title = "Red multicapa de transporte de Madrid (3D, enlaces entre capas)",
    scene = list(
      xaxis = list(title = "X (normalizada)"),
      yaxis = list(title = "Y (normalizada)"),
      zaxis = list(
        title = "Capa (modo de transporte)",
        tickvals = layer_index,
        ticktext = layer_order
      ),
      aspectmode = "cube"
    )
  )

  p3d

} else {
  message("Instala 'plotly'")
}
## Número de aristas inter-capa generadas: 728 
## Aristas con coordenadas válidas: 728 
## Puntos totales en edges_long: 2184

Interpretación de la visualización 3D multicapa

La visualización tridimensional obtenida permite observar de forma explícita la estructura multicapa del sistema de transporte urbano de Madrid, donde cada plano horizontal representa un modo de transporte (bicimad, bus, metro y parking) y las conexiones verticales indican zonas de intercambio multimodal.

Separación clara de capas y estructura espacial

En primer lugar, se aprecia una separación nítida entre las cuatro capas, lo que confirma que cada modo de transporte presenta una distribución espacial propia:

  • La capa bicimad aparece más concentrada espacialmente, reflejando su carácter urbano y centralizado.
  • La capa bus muestra una cobertura más extensa, coherente con su función de red capilar que conecta barrios periféricos.
  • La capa metro se sitúa entre ambas, con una estructura más compacta que refleja la jerarquía de estaciones y líneas.
  • La capa parking aparece claramente diferenciada y menos densa, lo que concuerda con su función de infraestructura de apoyo más localizada.

Esta separación valida la decisión metodológica de modelar el sistema como una red multicapa, ya que una red agregada ocultaría estas diferencias estructurales.

Uniones verticales como indicador de multimodalidad

Las líneas verticales que conectan nodos entre capas representan zonas donde coexisten varios modos de transporte. Estas conexiones no aparecen de forma homogénea, sino que se concentran en un subconjunto reducido de zonas, lo que pone de manifiesto una fuerte heterogeneidad en la multimodalidad del sistema.

En términos de teoría de redes:

  • Estas zonas actúan como nodos altamente versátiles, ya que participan simultáneamente en varias capas.

  • Son puntos críticos para la integración funcional del sistema, permitiendo el transbordo eficiente entre modos.

  • Su posición estructural las convierte en puntos de control del flujo global de movilidad.

Visualmente, estas zonas destacan como “columnas” que atraviesan varias capas, reforzando la idea de que la conectividad global del sistema depende de un número limitado de nodos clave.

Versatilidad y roles multinodales

La figura ilustra de manera intuitiva el concepto de versatilidad de nodo:

  • La mayoría de los nodos aparecen en una sola capa, actuando como nodos especializados (por ejemplo, paradas de bus o estaciones de bici aisladas).
  • En contraste, un número reducido de nodos conecta tres o incluso cuatro capas, desempeñando un rol multinodal esencial.

Este patrón es consistente con la literatura sobre redes multicapas, donde se observa que la importancia de un nodo no puede evaluarse correctamente desde una sola capa. La visualización 3D confirma que los nodos versátiles no solo existen, sino que son estructuralmente visibles cuando se representan explícitamente las capas y sus interdependencias.

Implicaciones para vulnerabilidad y resiliencia

Desde el punto de vista de la vulnerabilidad del sistema, la visualización sugiere una conclusión relevante:

  • La conectividad multimodal del transporte madrileño depende de un conjunto relativamente pequeño de nodos altamente interconectados entre capas.
  • Una perturbación localizada en estos nodos (por ejemplo, cierre de un intercambiador o fallo simultáneo de varios modos en una zona) podría propagarse rápidamente al resto del sistema, afectando a múltiples capas.

Este comportamiento es coherente con los modelos de cascadas de fallo en redes interdependientes descritos en la literatura (Buldyrev et al., 2010; Duan et al., 2019), donde la interdependencia aumenta tanto la eficiencia como la fragilidad del sistema.

Valor añadido de la representación multicapa

Finalmente, esta visualización demuestra el valor añadido del enfoque multicapa frente a representaciones tradicionales:

  • Permite identificar de forma directa los intercambiadores críticos.
  • Hace visible la estructura de interdependencia entre modos.
  • Facilita la interpretación de conceptos abstractos como versatilidad, roles multinodales y vulnerabilidad sistémica.

En conjunto, el gráfico 3D no solo actúa como herramienta exploratoria, sino como una evidencia visual clara de que el sistema de transporte de Madrid funciona como una red interconectada de capas, cuya eficiencia y resiliencia dependen de la correcta gestión de sus nodos multimodales clave.

5 Centralidad y versatilidad de las zonas

5.1 1. Medidas de centralidad por capa

Calculamos grado y fuerza (suma de pesos) para cada zona en cada modo.

centralidad_lista <- list()

for (modo in names(g_list)) {
  g <- g_list[[modo]]
  
  # Creamos el dataframe para la capa actual
  df_capa <- data.frame(
    cluster_id = as.integer(V(g)$name),
    modo       = modo,
    degree     = degree(g),
    strength   = strength(g, weights = E(g)$w),
    stringsAsFactors = FALSE
  )
  
  centralidad_lista[[modo]] <- df_capa
}

# Unimos todas las capas en un único dataframe
centralidad_por_capa <- do.call(rbind, centralidad_lista)

kable(head(centralidad_por_capa), caption = "Centralidad (grado y fuerza) por capa")
Centralidad (grado y fuerza) por capa
cluster_id modo degree strength
bicimad.1 1 bicimad 5 0.0071526
bicimad.2 2 bicimad 6 0.0131310
bicimad.3 3 bicimad 4 0.0072934
bicimad.4 4 bicimad 6 0.0105906
bicimad.5 5 bicimad 7 0.0130558
bicimad.6 6 bicimad 6 0.0151738

Distribución del grado por modo:

# Configuramos la rejilla de gráficos 2x2 para replicar el facet_wrap
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))

# Bucle para generar los histogramas por modo usando el vector 'modos'
for (m in modos) {
  # Filtrado manual sin pipes
  datos_sub <- centralidad_por_capa[centralidad_por_capa$modo == m, ]
  
  # Usamos el color asignado a cada modo en tu vector 'colores'
  hist(datos_sub$degree, 
       breaks = 30, 
       col = adjustcolor(colores[m], alpha.f = 0.6),
       border = "white",
       main = paste("MODO:", toupper(m)),
       xlab = "Grado",
       ylab = "Frecuencia")
}

# Restauramos el panel a 1x1
par(mfrow = c(1, 1))

5.2 2. Versatilidad de grado

En este trabajo utilizamos una medida de versatilidad de grado para caracterizar el papel multimodal de cada zona dentro de la red de transporte. Esta métrica no evalúa la intensidad del uso del sistema, sino la estructura de la conectividad de los nodos cuando se consideran simultáneamente las distintas capas modales (metro, autobús, cercanías, bicicleta pública, etc.).

La versatilidad de grado se define a partir de la entropía de Shannon de la distribución del grado de cada zona entre los distintos modos de transporte. Para cada nodo, el grado en cada capa se normaliza respecto al grado total, obteniendo una distribución de probabilidades que refleja cómo se reparte su conectividad entre modos. La entropía resultante se normaliza por su valor máximo teórico, de forma que la métrica toma valores entre 0 y 1.

Valores bajos de versatilidad indican zonas especializadas, cuya conectividad se concentra mayoritariamente en un único modo de transporte, incluso aunque su grado total pueda ser elevado. En cambio, valores altos corresponden a zonas multimodales, donde el grado se reparte de manera más equilibrada entre varias capas de la red. Estas zonas tienden a actuar como puntos de transferencia entre modos y desempeñan un papel clave en la articulación del sistema de transporte.

De este modo, la versatilidad de grado captura una propiedad estructural y relativa de la red, complementaria a las medidas tradicionales de centralidad. Mientras que el grado o el número de viajes reflejan la magnitud de la actividad, la versatilidad informa sobre la diversidad modal de dicha actividad y, por tanto, sobre la capacidad de una zona para redistribuir flujos entre distintos modos de transporte.

# --- 1. CREACIÓN DE LA MATRIZ DE GRADO ---
# Usamos xtabs para pivotar los datos de forma sencilla en R base
mat_grado_df <- as.data.frame.matrix(xtabs(degree ~ cluster_id + modo, data = centralidad_por_capa))

# Aseguramos que todas las columnas de modos existan y estén en orden
for (m in modos) {
  if (!m %in% names(mat_grado_df)) mat_grado_df[[m]] <- 0
}
grado_mat <- as.matrix(mat_grado_df[, modos])

# --- 2. CÁLCULO DE VERSATILIDAD ---
calc_versatilidad <- function(vec) {
  if (sum(vec) == 0) return(NA_real_)
  p <- vec / sum(vec)
  p <- p[p > 0]
  H <- -sum(p * log(p))
  # Normalización por el logaritmo del número total de capas
  H / log(length(vec))
}

versatilidad <- apply(grado_mat, 1, calc_versatilidad)

# --- 3. PREPARACIÓN DE UNIONES (JOINS) ---
# Aseguramos tipos de datos para las claves de unión
zonas_attrs$distrito_asignado <- as.integer(zonas_attrs$distrito_asignado)
viajes_distrito$distrito <- as.integer(viajes_distrito$distrito)

# Creamos el dataframe base de versatilidad
versatilidad_df <- data.frame(
  cluster_id = as.integer(rownames(grado_mat)),
  versatilidad_grado = as.numeric(versatilidad),
  stringsAsFactors = FALSE
)

# Realizamos los left_join usando merge
versatilidad_df <- merge(versatilidad_df, zonas_attrs, by = "cluster_id", all.x = TRUE)
versatilidad_df <- merge(versatilidad_df, viajes_distrito, 
                         by.x = "distrito_asignado", 
                         by.y = "distrito", 
                         all.x = TRUE)

# --- 4. RESULTADOS ---
print("Columnas en versatilidad_df:")
## [1] "Columnas en versatilidad_df:"
print(colnames(versatilidad_df))
##  [1] "distrito_asignado"         "cluster_id"               
##  [3] "versatilidad_grado"        "stop_lat"                 
##  [5] "stop_lon"                  "transport_mode"           
##  [7] "bicimad"                   "bus"                      
##  [9] "metro"                     "parking"                  
## [11] "personas"                  "viajes_totales"           
## [13] "viajes_medios_por_persona" "personas_totales"
summary(versatilidad_df$versatilidad_grado)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.1477  0.0000  0.9983

5.2.1 Nodos más versátiles y más especializados

# --- 1. TOP 15 ZONAS MÁS VERSÁTILES ---
orden_desc <- order(versatilidad_df$versatilidad_grado, decreasing = TRUE)
top_versatiles_full <- head(versatilidad_df[orden_desc, ], n = 15)

# Verificamos qué columnas están disponibles para evitar el error "undefined columns"
cols_deseadas <- c("cluster_id", "versatilidad_grado", "modos_disponibles", "distrito_asignado")
cols_presentes <- cols_deseadas[cols_deseadas %in% colnames(top_versatiles_full)]

top_versatiles <- top_versatiles_full[, cols_presentes]

# --- 2. TOP 15 ZONAS MÁS ESPECIALIZADAS ---
orden_asc <- order(versatilidad_df$versatilidad_grado, na.last = TRUE)
top_especializados_full <- head(versatilidad_df[orden_asc, ], n = 15)

top_especializados <- top_especializados_full[, cols_presentes]

# --- 3. TABLAS RESULTADO ---
kable(top_versatiles, caption = "Top 15 zonas más versátiles")
Top 15 zonas más versátiles
cluster_id versatilidad_grado distrito_asignado
265 164 0.9983119 2807907
108 446 0.9976183 2807904
286 101 0.9968869 2807907
264 41 0.9955380 2807907
96 283 0.9949672 2807903
101 56 0.9949672 2807904
5 239 0.9946468 2807901
112 814 0.9934804 2807904
117 71 0.9927376 2807904
407 534 0.9927376 2807908
738 339 0.9927376 2807911
29 7 0.9926388 2807901
1173 251 0.9892423 2807916
1432 562 0.9892423 2807920
124 573 0.9863237 2807904
kable(top_especializados, caption = "Top 15 zonas más especializadas")
Top 15 zonas más especializadas
cluster_id versatilidad_grado distrito_asignado
1 929 0 2807901
2 732 0 2807901
4 1410 0 2807901
10 1500 0 2807901
11 674 0 2807901
12 731 0 2807901
17 1503 0 2807901
23 140 0 2807901
25 1488 0 2807901
31 530 0 2807901
32 890 0 2807901
33 673 0 2807901
34 902 0 2807901
37 1526 0 2807901
43 1135 0 2807902

Desde el punto de vista espacial, las zonas con alta versatilidad de grado se concentran en torno a los grandes intercambiadores de transporte y a las áreas más centrales de la ciudad. Se trata de nodos que combinan conectividad en varias capas de la red (metro, bus, cercanías, bicicleta pública), de modo que actúan como puntos de transferencia entre modos diferentes. Estas zonas tienen grados relativamente altos en varias capas simultáneamente y, por tanto, desempeñan un papel de “puente multimodal” entre distintas partes del sistema de transporte.

En el extremo opuesto, las zonas con baja versatilidad tienden a estar especializadas en un único modo de transporte, típicamente la red de autobuses o una línea concreta de metro. Su grado puede ser elevado dentro de una sola capa, pero su conectividad se limita a ese modo, de modo que su capacidad para redistribuir flujos entre capas es reducida. Esta especialización puede ser eficiente en términos operativos, pero hace que la zona sea menos flexible ante perturbaciones que afecten a su modo dominante, ya que existen menos alternativas multimodales inmediatas.

5.3 3. Relación entre versatilidad e intensidad de viajes

# Quitamos NAs de las columnas que vamos a usar
datos_plot <- na.omit(versatilidad_df[, c("viajes_medios_por_persona", "versatilidad_grado")])

# Ordenamos por el eje X (necesario para que la línea no se cruce)
datos_plot <- datos_plot[order(datos_plot$viajes_medios_por_persona), ]

# Gráfico base
plot(datos_plot$viajes_medios_por_persona, datos_plot$versatilidad_grado,
     pch = 19, col = adjustcolor("black", 0.3),
     xlab = "Viajes medios por persona", ylab = "Versatilidad")

# Línea de tendencia (loess)
modelo <- loess(versatilidad_grado ~ viajes_medios_por_persona, data = datos_plot)
lines(datos_plot$viajes_medios_por_persona, predict(modelo), col = "blue", lwd = 2)

La Figura anterior muestra la relación entre la y el a nivel de distrito. Visualmente, la nube de puntos no presenta una tendencia clara, lo que sugiere una relación débil entre ambas variables. Este resultado se confirma mediante el análisis de correlación de Pearson.

datos_cor <- versatilidad_df %>% 
  filter(!is.na(versatilidad_grado),
         !is.na(viajes_medios_por_persona))

n_obs <- nrow(datos_cor)
n_obs
## [1] 1558
if (n_obs >= 3) {
  cor_test <- cor.test(
    datos_cor$versatilidad_grado,
    datos_cor$viajes_medios_por_persona
  )
  cor_test
} else {
  cat("No hay suficientes observaciones completas para estimar de forma fiable la correlación entre versatilidad y viajes.
")
}
## 
##  Pearson's product-moment correlation
## 
## data:  datos_cor$versatilidad_grado and datos_cor$viajes_medios_por_persona
## t = 0.14253, df = 1556, p-value = 0.8867
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.04605711  0.05326588
## sample estimates:
##         cor 
## 0.003613297

Este cálculo arroja un coeficiente prácticamente nulo (\(r = 0.0036\)) y no significativo (\(p = 0.8867\), \(n = 1558\)). El intervalo de confianza al 95,% incluye el valor cero, lo que indica la ausencia de una asociación lineal estadísticamente significativa entre ambas magnitudes.

La línea de suavizado refuerza esta interpretación, mostrando un comportamiento prácticamente plano a lo largo del rango de valores observados. Aunque se aprecia una ligera curvatura en los valores más altos de viajes medios, esta variación es reducida y carece de relevancia estadística. En consecuencia, no se observa evidencia de que los distritos con mayor intensidad de movilidad presenten zonas con una conectividad más equilibrada entre modos de transporte.

Este resultado es coherente con la definición de la versatilidad de grado como una medida basada en la distribución del grado entre capas modales, y no en el volumen absoluto de actividad. Una zona puede registrar un número elevado de viajes concentrados en un único modo y, por tanto, mostrar una baja versatilidad, mientras que otras zonas con menor intensidad de uso pueden desempeñar un papel multimodal más destacado. De este modo, la versatilidad de grado captura una dimensión complementaria de la red, asociada a su diversidad modal, que no depende directamente de la cantidad de viajes generados.

6 Vulnerabilidad de la red de transporte

En lugar de implementar un modelo muy complejo de cascadas, realizamos un análisis sencillo pero informativo de robustez estructural:

  • Tomamos la red agregada (unión de todas las capas).
  • Eliminamos nodos de forma dirigida (los más centrales) en una capa y observamos cómo disminuye el tamaño de la componente gigante agregada.

6.1 1. Red agregada

# Matrices de adyacencia por capa
adj_list <- lapply(g_list, function(g) {
  as.matrix(as_adj(g, attr = "w", sparse = FALSE))
})

# Aseguramos mismo orden de nodos
orden_nodos <- order(as.integer(V(g_list[[1]])$name))
adj_list <- lapply(adj_list, function(m) m[orden_nodos, orden_nodos])

adj_agregada <- Reduce("+", adj_list)

g_agregado <- graph_from_adjacency_matrix(adj_agregada, mode = "undirected", weighted = TRUE)

comp <- components(g_agregado)
tamano_gcc_inicial <- max(comp$csize)
tamano_gcc_inicial
## [1] 1528

6.2 2. Robustez frente a fallos en una capa

Analizamos qué ocurre si vamos eliminando nodos con mayor grado en la capa metro y miramos el tamaño relativo de la componente gigante en la red agregada.

# --- 1. PREPARACIÓN DEL ORDEN DE FALLO ---
g_metro <- g_list[["metro"]]
deg_metro <- degree(g_metro)

# Nombres de los nodos ordenados de mayor a menor grado
orden_fallo <- names(sort(deg_metro, decreasing = TRUE))

# --- 2. SIMULACIÓN DE ROBUSTEZ ---
fracciones <- seq(0, 0.5, by = 0.05)
res_lista <- list()

for (i in seq_along(fracciones)) {
  f <- fracciones[i]
  
  # Calculamos cuántos nodos eliminar
  n_remove <- ceiling(f * length(orden_fallo))
  nodos_fuera <- if (n_remove > 0) orden_fallo[1:n_remove] else character(0)
  
  # Eliminamos los nodos de la red agregada y calculamos el GCC
  g_agregado_reducido <- delete_vertices(g_agregado, nodos_fuera)
  comp_red <- components(g_agregado_reducido)
  
  gcc_red <- if (length(comp_red$csize) == 0) 0 else max(comp_red$csize)
  
  # Guardamos resultados en un dataframe temporal
  res_lista[[i]] <- data.frame(
    frac_eliminada = f,
    tamano_gcc_rel = gcc_red / tamano_gcc_inicial
  )
}

# Unimos los resultados
resultado_robustez <- do.call(rbind, res_lista)

# --- 3. VISUALIZACIÓN ---
# Configuramos el gráfico base
plot(resultado_robustez$frac_eliminada, resultado_robustez$tamano_gcc_rel,
     type = "b", # Tipo "both" (puntos y líneas)
     pch = 19, 
     col = "#2c7bb6", # Azul Metro
     lwd = 2,
     main = "Robustez ante fallos en la capa METRO",
     xlab = "Fracción de nodos eliminados (%)",
     ylab = "Tamaño relativo del GCC (%)",
     xaxt = "n", yaxt = "n",
     frame.plot = FALSE)

# Personalizamos los ejes para mostrar porcentajes
eje_x <- seq(0, 0.5, 0.1)
axis(1, at = eje_x, labels = paste0(eje_x * 100, "%"))

eje_y <- seq(0, 1, 0.2)
axis(2, at = eje_y, labels = paste0(eje_y * 100, "%"), las = 1)

grid(col = "gray90", lty = "solid")

La Figura muestra la evolución del tamaño relativo de la componente gigante de la red agregada a medida que se eliminan progresivamente nodos de la capa de metro, siguiendo un criterio dirigido basado en su centralidad. Se observa una disminución aproximadamente monótona y casi lineal de la conectividad global, lo que indica que los nodos más centrales de la red de metro desempeñan un papel estructural clave en la cohesión del sistema multicapa.

En particular, la eliminación del 20–30 % de los nodos más relevantes de la capa de metro reduce el tamaño de la componente gigante a alrededor del 70–75 % de su valor inicial, mientras que al eliminar la mitad de los nodos centrales la red pierde más del 50 % de su conectividad global. Este resultado pone de manifiesto que la capa de metro actúa como columna vertebral del sistema de transporte urbano: disrupciones severas concentradas en sus nodos más importantes pueden propagarse al resto de capas y provocar una fragmentación significativa de la red agregada.

6.2.1 Comparación entre fallos dirigidos y aleatorios

set.seed(123)

simular_ataque <- function(g, estrategia = c("dirigido", "aleatorio"),
                           fracciones = seq(0, 0.5, by = 0.05)) {
  estrategia <- match.arg(estrategia)
  n <- vcount(g)
  resultados <- tibble()
  
  if (estrategia == "dirigido") {
    orden <- names(sort(degree(g), decreasing = TRUE))
  } else {
    orden <- sample(V(g)$name)
  }
  
  for (f in fracciones) {
    k <- ceiling(f * n)
    g_tmp <- delete_vertices(g, orden[seq_len(k)])
    gcc <- if (vcount(g_tmp) == 0) 0 else max(components(g_tmp)$csize)
    
    resultados <- bind_rows(
      resultados,
      tibble(
        frac_eliminada = f,
        tamano_gcc_rel = gcc / tamano_gcc_inicial,
        estrategia = estrategia
      )
    )
  }
  resultados
}

res_dirigido  <- simular_ataque(g_agregado, "dirigido")
res_aleatorio <- simular_ataque(g_agregado, "aleatorio")

res_vuln <- bind_rows(res_dirigido, res_aleatorio)

ggplot(res_vuln,
       aes(x = frac_eliminada, y = tamano_gcc_rel, color = estrategia)) +
  geom_line(size = 1) +
  geom_point() +
  theme_minimal() +
  labs(
    x = "Fracción de nodos eliminados",
    y = "Tamaño relativo de la componente gigante",
    color = "Estrategia",
    title = "Robustez de la red: ataque dirigido vs fallo aleatorio"
  )

El análisis conjunto de la versatilidad nodal y la robustez estructural de la red agregada de Madrid arroja tres conclusiones fundamentales sobre la arquitectura del sistema de transporte:

  1. Independencia entre Infraestructura y Demanda Individual. Tal como se observa en el gráfico de Versatilidad vs. Viajes, no existe una correlación clara entre la intensidad de uso (viajes medios por persona en el distrito) y la complejidad multimodal de las paradas.

    • Estratificación de la oferta: La disposición de los puntos en bandas horizontales discretas revela que la red no es continua, sino jerárquica: la gran mayoría de nodos son unimodales (base del gráfico), mientras que la multimodalidad se da “a saltos” (agregando capas de Metro o BiciMAD).

    • Homogeneidad de la demanda: El rango estrecho del eje X (1.95 - 2.15 viajes) confirma que la movilidad individual es uniforme en toda la ciudad. Por tanto, la ubicación de los grandes intercambiadores (nodos con alta versatilidad) responde a criterios estratégicos de planificación de red y transbordo, y no a una mayor intensidad de viajes per cápita en esos distritos específicos.

  2. La red es “Robusta pero Frágil”. La comparación entre estrategias de fallo (gráfico de Ataque dirigido vs. Fallo aleatorio) confirma que la red de transporte de Madrid comparte la topología típica de las redes complejas libres de escala (scale-free):

    • Alta tolerancia a fallos aleatorios (Línea Roja): El sistema es extremadamente resiliente ante incidencias fortuitas. Incluso ante la inoperatividad del 40% de los nodos de forma aleatoria (simulando averías dispersas o cortes menores), la red mantiene una componente gigante conectada significativa, garantizando la movilidad en gran parte del territorio.

    • Extrema fragilidad ante ataques dirigidos (Línea Azul/Turquesa): Por el contrario, la red es muy vulnerable si los fallos afectan a sus nodos centrales. La eliminación de apenas el 15-20% de los nodos más conectados provoca el colapso funcional del sistema (fragmentación de la componente gigante por debajo del 20%), aislando grandes zonas de la ciudad.

  3. El Metro como columna vertebral crítica. El análisis específico de la capa de Metro (gráfico de caída monótona) demuestra su rol de esqueleto estructural. Al eliminar los nodos principales de esta capa, la conectividad global de la red agregada decae casi linealmente. Esto implica que la cohesión entre los distintos modos (bus, bici, parking) depende críticamente de la operatividad de las estaciones de metro centrales, que actúan como los “pegamentos” que unen las distintas capas de la red multicapa.

7 Dashboard Interactivo de la Red de Transporte

En esta sección se presenta un dashboard interactivo que resume las métricas principales del análisis de redes multicapas de Madrid. Este dashboard permite una exploración visual e interactiva de los resultados obtenidos.

7.1 Métricas Generales del Sistema

7.2 Distribución de Nodos por Capa

7.3 Top 10 Zonas Más Versátiles

7.4 Análisis de Centralidad por Capa

7.5 Robustez: Evolución de la Componente Gigante

7.6 Resumen Ejecutivo

Conclusiones Clave del Dashboard

  • Sistema Integrado: La red de Madrid cuenta con 1558 zonas de intercambio que integran 3573 conexiones entre diferentes modos de transporte.
  • Predominancia del Bus: La capa de bus representa el modo de transporte con mayor número de nodos, seguido por el metro, reflejando la capilaridad del sistema.
  • Zonas Estratégicas: Las zonas con mayor versatilidad (top 10) actúan como hubs críticos que conectan múltiples capas de transporte.
  • Vulnerabilidad Diferencial: La red muestra alta robustez ante fallos aleatorios pero extrema fragilidad ante ataques dirigidos a nodos centrales, especialmente en la capa de metro.
  • Movilidad Homogénea: Los viajes medios por persona son notablemente uniformes entre distritos (≈2 viajes/día), independientemente de la versatilidad de las zonas.

8 Conclusiones

En esta memoria se ha construido y analizado una red multicapas de transporte para Madrid utilizando exclusivamente los datos proporcionados:

  • Las paradas se han agrupado en zonas de intercambio mediante clustering espacial.
  • Se ha definido una red multicapa con capas para bus, metro, bicimad y parking.
  • Se ha calculado la centralidad de las zonas en cada capa y un indicador de versatilidad de grado que muestra qué zonas son importantes en varios modos a la vez.
  • Se ha generado una visualización multicapa con muxViz que hace explícita la estructura en capas.
  • Se ha estudiado la robustez estructural de la red agregada frente a fallos dirigidos en nodos centrales de la capa metro.

A partir de aquí, se podrían ampliar los análisis:

  • Incorporando información de líneas reales (no solo proximidad espacial).
  • Usando otras medidas de centralidad (betweenness, eigenvector, PageRank).
  • Implementando modelos de cascadas más cercanos a la literatura de redes interdependientes.

El código está preparado para ejecutarse directamente (asumiendo que los CSV y los paquetes necesarios están disponibles) y centrado únicamente en los datos reales de Madrid, tal y como se pedía en el enunciado.

8.1 Implementación Técnica y Uso de muxViz

Para la realización de este estudio, la librería muxViz en R ha sido la herramienta fundamental, permitiéndonos traducir los datos crudos de transporte en un modelo matemático robusto. El proceso se ha estructurado en tres fases clave:

Modelado de la Red: A partir de las listas de enlaces procesadas, utilizamos la función buildMultilayerNetworkFromEdgeList(). Esto fue crucial para generar el objeto multicapa (mlo) y construir internamente la supra-matriz de adyacencia, que unifica las capas de Metro, Bus y BiciMAD en una sola estructura tensorial.

Análisis de Métricas: Para cuantificar la relevancia de las paradas, calculamos el Grado Multicapa mediante GetMultiLayerDegree(), obteniendo la conectividad total. Sin embargo, para medir la versatilidad real de los nodos (su capacidad para influir en todo el sistema y no solo en su capa), nos basamos en el cálculo del MultiLayer PageRank (GetMultiLayerPageRank()), identificando así los verdaderos hubs intermodales de Madrid.

Visualización: Finalmente, para interpretar la topología resultante, empleamos layout_multiplex() para fijar las coordenadas geoespaciales y plot_multiplex() para generar las representaciones gráficas en 3D, permitiendo observar visualmente la interacción y el solapamiento entre las distintas capas de transporte.